/**
 * @author Ahammer
 * @date   2009 07
 * @update 2010 01 added progress bar support
 * @update 2010 03 added grey value image support, Surface is calculated according to Chinga etal. J Microsc. 227 2007
 *                 not sure if it works properly!!!
 * @update 2011 01 bug fix for binary images. Condition for filteredsubsample : resizing previously resized image	
 * @update 2011 01 added 16 bit support 
 * @update 2011 02 added regression data to table 
 * @update 2011 02 removing ModelListener extremely speeds up model configuration (adding columns, writing data,...)                 
 */

/*
 * This file is part of Iqm.
 * Copyright (c) 2010-2011 Helmut Ahammer
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package op;

import java.awt.RenderingHints;
import java.awt.image.DataBuffer;
import java.awt.image.Raster;
import java.awt.image.renderable.ParameterBlock;
import java.util.Vector;
import javax.media.jai.BorderExtender;
import javax.media.jai.Histogram;
import javax.media.jai.Interpolation;
import javax.media.jai.JAI;
import javax.media.jai.KernelJAI;
import javax.media.jai.ParameterBlockJAI;
import javax.media.jai.PlanarImage;
import javax.media.jai.RenderedOp;
import javax.swing.JTable;
import javax.swing.table.DefaultTableModel;
import main.IqmTools;
import main.Manager;

/**
 * This is the main image processing class
 * There exist two approaches:
 * A user defined JAI operator is just called 
 * or:
 * The actual processing is implemented in this class 
 */
public class IqmFracPyramidOperator{ 


  public IqmFracPyramidOperator() {
	  //WARNING: Don't declare fields here
	  //Fields declared here aren't thread safe!   
  }
  /**
   * 
   * @param ParametrBlockJAI pb
   */
  public Object run(ParameterBlockJAI pbJAI){
	  	Object ob = null;
	    //Object ob = JAI.create("IqmFracPyramid".toLowerCase(), pb, null); //does not work
	    //therefore, the implementation is here:
	  	PlanarImage pi = (PlanarImage) pbJAI.getSource(0);
	  	int number     = pbJAI.getIntParameter("NumPyrImgs");
	  	int optIntP    = pbJAI.getIntParameter("Interpolation");
	  	int regStart   = pbJAI.getIntParameter("RegStart");
	  	int regEnd     = pbJAI.getIntParameter("RegEnd");
	  	boolean optShowPlot = false;
	  	boolean optDeleteExistingPlot = false;
	  	if (pbJAI.getIntParameter("ShowPlot") == 1 )          optShowPlot           = true;	  	
	  	if (pbJAI.getIntParameter("DeleteExistingPlot") == 1) optDeleteExistingPlot = true;	  	

	  	int numBands  = pi.getData().getNumBands();
	  	String type = IqmTools.getImgTyp(pi);
	    double typeGreyMax = IqmTools.getImgTypeGreyMax(pi);  	
	  	boolean isBinary = true;
//	  	int pixelSize = pi.getColorModel().getPixelSize();
//		String type = IqmTools.getCurrImgTyp(pixelSize, numBands);
	  	String imgName  = (String) pi.getProperty("image_name");
	  	String fileName = (String) pi.getProperty("file_name");
 
	  	//initialize table for Dimension data
		DefaultTableModel model = new DefaultTableModel();
		JTable jTable  = new JTable(model);	
		//adding a lot of columns would be very slow due to active model listener
		model.removeTableModelListener(jTable);
		model.addColumn("FileName");
		model.addColumn("ImageName"); 
		model.addColumn("Band");
		model.addColumn("RegStart"); 
		model.addColumn("RegEnd"); 
		model.addColumn("Interp"); 
		for(int b = 0; b < pi.getNumBands(); b++){ //mehrere Bands
			model.addRow(new Object[]{fileName, imgName, b, regStart, regEnd, optIntP});
		}	
	  	double[][] totals    = new double[number][numBands];
	  	double[]   totalsMax = new double[numBands];
	  	double[][] eps       = new double[number][numBands];
	  	
	  	//definition of eps
		for (int n = 0; n < number; n++){
			for(int b = 0; b < numBands; b++){
				eps[n][b] = 1.0/Math.pow(2, n);
				System.out.println("IqmFracPyramidOperator: n:" + n + " eps:  "+  eps[n][b] );
			}
		}			
	  	
	  	//first  image
	    // Set up the parameters for the Histogram object.
	    int[]    bins  = {(int)typeGreyMax+1, (int)typeGreyMax+1, (int)typeGreyMax+1}; // The number of bins  e.g. {256, 256, 256}
	    double[] lows  = {0.0D, 0.0D, 0.0D};                                           // The low incl.value  e.g. {0.0D, 0.0D, 0.0D}
	    double[] highs = {typeGreyMax+1, typeGreyMax+1, typeGreyMax+1};                // The high excl.value e.g. {256.0D, 256.0D, 256.0D}
	    // Create the parameter block.
	    ParameterBlock pbHisto = new ParameterBlock();
	    pbHisto.addSource(pi);   //Source image
	    pbHisto.add(null);       //ROI
	    pbHisto.add(1);          //Horizontal sampling
	    pbHisto.add(1);          //Vertical sampling
	    pbHisto.add(bins);       //Number of Bins
	    pbHisto.add(lows);       //Lowest inclusive values
	    pbHisto.add(highs);      //Highest exclusive values
	    // Perform the histogram operation.
		RenderedOp rOp = JAI.create("Histogram", pbHisto, null);
		// Retrieve the histogram data.
		Histogram histogram = (Histogram)rOp.getProperty("histogram");
		//System.out.println("IqmFracPyramidOperator: (int)typeGreyMax: " + (int)typeGreyMax);
		double totalGrey   = (double)histogram.getSubTotal(0, 1, (int)typeGreyMax);  //without 0!
		double totalBinary = (double)histogram.getSubTotal(0, (int)typeGreyMax, (int)typeGreyMax);
		if (totalBinary == totalGrey){
			isBinary = true;
			System.out.println("IqmFracPyramidOperator: Binary image = true");
		} else{
			isBinary = false;
			System.out.println("IqmFracPyramidOperator: Binary image = false");
		}	
		for (int b = 0; b < numBands; b++) {
			//totals[0][b] = (double)histo.getSubTotal(b, 1, (int)typeGreyMax);  //without 0!
			if (isBinary) {
				totalBinary = (double)histogram.getSubTotal(b, (int)typeGreyMax, (int)typeGreyMax); //[1, (int)typeGreyMax] or [(int)typeGreyMax, (int)typeGreyMax] both ok
				totalsMax[b] = totalBinary; //Binary Image
				totals[0][b] = totalBinary/totalsMax[b]; //Binary Image
			}else{
				totalsMax[b] = pi.getWidth() * pi.getHeight(); 
				//Surface according to Chinga etal. J Microsc. 227 2007
				ParameterBlock pbF = new ParameterBlock();
			    pbF.addSource(pi);
				pbF.add(DataBuffer.TYPE_FLOAT);
				PlanarImage piF = JAI.create("Format", pbF, null);
		
				float[] h_data = { 1.0f,  0.0f, -1.0f,
					  			   2.0f,  0.0f, -2.0f,
					  			   1.0f,  0.0f, -1.0f  };
				float[] v_data = {-1.0f, -2.0f, -1.0f,
							       0.0f,  0.0f,  0.0f,
							       1.0f,  2.0f,  1.0f  };
				for (int i = 0; i < h_data.length; i++){
					h_data[i] = h_data[i]*(1.0f/8.0f);
				}
				for (int i = 0; i < v_data.length; i++){
					v_data[i] = v_data[i]*(1.0f/8.0f);
				}
				KernelJAI kernel_h = new KernelJAI(3,3, h_data);
				KernelJAI kernel_v = new KernelJAI(3,3, v_data);
				RenderingHints rh = new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY));	
							
				ParameterBlock pbGH = new ParameterBlock();
			    pbGH.addSource(piF);
				pbGH.add(kernel_h);
				PlanarImage piGH = JAI.create("Convolve", pbGH, rh);
				
				ParameterBlock pbM = new ParameterBlock();
				pbM.addSource(piGH);
				pbM.addSource(piGH);
				PlanarImage piGH2 = JAI.create("Multiply", pbM, null);
				
				ParameterBlock pbGV = new ParameterBlock();
			    pbGV.addSource(piF);
				pbGV.add(kernel_v);
				PlanarImage piGV = JAI.create("Convolve", pbGV, rh);
				
				pbM = new ParameterBlock();
				pbM.addSource(piGV);
				pbM.addSource(piGV);
				PlanarImage piGV2 = JAI.create("Multiply", pbM, null);
				
				ParameterBlock pbA = new ParameterBlock();
				pbA.addSource(piGH2);
				pbA.addSource(piGV2);
				PlanarImage piGM = JAI.create("Add", pbA, null);
			
				ParameterBlock pbGM = new ParameterBlock();
				pbGM.addSource(piGM);
				pbGM.add(new double[]{1.0d, 1.0d, 1.0d});
				piGM = JAI.create("AddConst", pbGM, null);
				
				//calculate square root
				ParameterBlock pbLog = new ParameterBlock();
				pbLog.addSource(piGM);
				PlanarImage piLog = JAI.create("Log", pbLog, null);
				
				ParameterBlock pbMConst = new ParameterBlock();
				pbMConst.addSource(piLog);
				pbMConst.add(new double[]{0.5d, 0.5d, 0.5d});
				piLog = JAI.create("MultiplyConst", pbMConst, null);
				
				ParameterBlock pbEx = new ParameterBlock();
				pbEx.addSource(piLog);
				PlanarImage piEx = JAI.create("Exp", pbEx, null);
								
				Raster r = piEx.getData();
				int minX = r.getMinX();
				int minY = r.getMinY();
				for(int x = 0; x < r.getWidth(); x++){
					for(int y = 0; y < r.getHeight(); y++){
						totals[0][b]= totals[0][b] + r.getSample(minX + x, minY + y, b) ;
					}
				}
				//totals[0][b] = totals[0][b] / totalsMax[b];
			} //if isbinary == false
			
		} // b band		 
	  		
		PlanarImage pi2 = pi; //pi2 will be the size reduced image

		for (int n = 1; n < number; n++){
			int proz = (n+1)*100 / number;
	   		Manager.setProgressBar(proz);	
			if ((optIntP == 0) || (optIntP == 1) || (optIntP == 2) ||(optIntP == 3) ){
				//condition for filteredsubsample : resize previously resized image	
				ParameterBlock pb = new ParameterBlock();
				pb.addSource(pi2);
				pb.add(2);//x downsample factor
				pb.add(2);//y downsample factor
				//pb.add((int)Math.pow(2, n));//x downsample factor
				//pb.add((int)Math.pow(2, n));//y downsample factor
				
				//Gauss filter
				pb.add(new float[] {1.0f}); //for binary images 1.0f
				//pb.add(new float[] {0.0632f, 0.0558f});
				//pb.add(new float[] {0.0632f, 0.0558f, 0.0383f});
				//pb.add(new float[] {0.0117f, 0.0115f, 0.0108f, 0.0098f, 0.0085f, 0.0071f});
				
				if (optIntP == 0) pb.add(Interpolation.getInstance(Interpolation.INTERP_NEAREST));
				if (optIntP == 1) pb.add(Interpolation.getInstance(Interpolation.INTERP_BILINEAR));
				if (optIntP == 2) pb.add(Interpolation.getInstance(Interpolation.INTERP_BICUBIC));
				if (optIntP == 3) pb.add(Interpolation.getInstance(Interpolation.INTERP_BICUBIC_2));
			
				RenderingHints rh = new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY));
				pi2 = JAI.create("filteredsubsample", pb, rh);
				//pi2 = JAI.create("scale", pb, rh);
				
			}
			if (optIntP == 4) {	
				//condition for subsampleaverage : resize always from original sized image!
				ParameterBlock pb = new ParameterBlock();
				pb.addSource(pi);
				pb.add(1.0/Math.pow(2, n));//x downsample factor
				pb.add(1.0/Math.pow(2, n));//y downsample factor


				RenderingHints rh = new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY));
				//rh = null;
				pi2 = JAI.create("subsampleaverage", pb, rh);
				//System.out.println("IqmFracPyramidOperator: pi2.getWidth "+ pi2.getWidth());
			}
			//create binary image, not necessary
//			double[] low  = new double[numBands];
//			double[] high = new double[numBands];
//			double[] map  = new double[numBands];			
//			for (int i = 0; i < numBands; i++){
//				low[i]  = 1d;
//				high[i] = (int)typeGreyMax-1;
//				map[i]  = (int)typeGreyMax;
//			}
//			ParameterBlock pbBinary = new ParameterBlock();
//			pbBinary.addSource(pi2);
//			pbBinary.add(low);
//			pbBinary.add(high);
//			pbBinary.add(map);		 
//			pi2 = JAI.create("threshold", pbBinary);
					
			//Optional wird eine Bild ausgegeben
//			if (n == 3){
//				ob = pi2;
//			}
			
			if (isBinary) {
				pbHisto.removeSources();
			    pbHisto.addSource(pi2);   //Source image
			    // Perform the histogram operation.
				rOp = JAI.create("Histogram", pbHisto, null);
				// Retrieve the histogram data.
				histogram = (Histogram)rOp.getProperty("histogram");
			
			}else{
				//rOp = JAI.create("Extrema", pi2, null);			
			}
	
			for (int b = 0; b < numBands; b++) {
				//totals[n][b] = (double) histo.getSubTotal(b, 1, (int)typeGreyMax-1); only binary
				if (isBinary) {
					totals[n][b] = (double) histogram.getSubTotal(b, 1, (int)typeGreyMax); //Binary Image //[1, 255] and not [255, 255] because interpolation introduces grey values other than 255!
					totals[n][b] = totals[n][b] / totalsMax[b];	
				}else{
					//Surface according to Chinga etal. J Microsc. 227 2007
					ParameterBlock pbF = new ParameterBlock();
				    pbF.addSource(pi2);
					pbF.add(DataBuffer.TYPE_FLOAT);
					PlanarImage piF = JAI.create("Format", pbF, null);
			
					float[] h_data = { 1.0f,  0.0f, -1.0f,
						  			   2.0f,  0.0f, -2.0f,
						  			   1.0f,  0.0f, -1.0f  };
					float[] v_data = {-1.0f, -2.0f, -1.0f,
								       0.0f,  0.0f,  0.0f,
								       1.0f,  2.0f,  1.0f  };
					for (int i = 0; i < h_data.length; i++){
						h_data[i] = h_data[i]*(1.0f/8.0f);
					}
					for (int i = 0; i < v_data.length; i++){
						v_data[i] = v_data[i]*(1.0f/8.0f);
					}
					KernelJAI kernel_h = new KernelJAI(3,3, h_data);
					KernelJAI kernel_v = new KernelJAI(3,3, v_data);
					RenderingHints rh = new RenderingHints(JAI.KEY_BORDER_EXTENDER, BorderExtender.createInstance(BorderExtender.BORDER_COPY));	
								
					ParameterBlock pbGH = new ParameterBlock();
				    pbGH.addSource(piF);
					pbGH.add(kernel_h);
					PlanarImage piGH = JAI.create("Convolve", pbGH, rh);
					
					ParameterBlock pbM = new ParameterBlock();
					pbM.addSource(piGH);
					pbM.addSource(piGH);
					PlanarImage piGH2 = JAI.create("Multiply", pbM, null);
					
					ParameterBlock pbGV = new ParameterBlock();
				    pbGV.addSource(piF);
					pbGV.add(kernel_v);
					PlanarImage piGV = JAI.create("Convolve", pbGV, rh);
					
					pbM = new ParameterBlock();
					pbM.addSource(piGV);
					pbM.addSource(piGV);
					PlanarImage piGV2 = JAI.create("Multiply", pbM, null);
					
					ParameterBlock pbA = new ParameterBlock();
					pbA.addSource(piGH2);
					pbA.addSource(piGV2);
					PlanarImage piGM = JAI.create("Add", pbA, null);
				
					ParameterBlock pbGM = new ParameterBlock();
					pbGM.addSource(piGM);
					pbGM.add(new double[]{1.0d, 1.0d, 1.0d});
					piGM = JAI.create("AddConst", pbGM, null);
					
					//Square root
					ParameterBlock pbLog = new ParameterBlock();
					pbLog.addSource(piGM);
					PlanarImage piLog = JAI.create("Log", pbLog, null);
					
					ParameterBlock pbMConst = new ParameterBlock();
					pbMConst.addSource(piLog);
					pbMConst.add(new double[]{0.5d, 0.5d, 0.5d});
					piLog = JAI.create("MultiplyConst", pbMConst, null);
					
					ParameterBlock pbEx = new ParameterBlock();
					pbEx.addSource(piLog);
					PlanarImage piEx = JAI.create("Exp", pbEx, null);
					
					
//					rOp = JAI.create("Histogram", piGM);			
//					histo = (Histogram)rOp.getProperty("histogram");
//					totals[n][b] = histo.getTotals()[b];
					
//					ParameterBlock pbGM = new ParameterBlock();
//				    pbGM.addSource(piF);
//					pbGM.add(kernel_h);
//					pbGM.add(kernel_v);
//					PlanarImage piGM = JAI.create("GradientMagnitude", pbGM, rh);
					
					Raster r = piEx.getData();
					int minX = r.getMinX();
					int minY = r.getMinY();
					for(int x = 0; x < r.getWidth(); x++){
						for(int y = 0; y < r.getHeight(); y++){
							//System.out.println("IqmFracOperator:  r.getSample(minX + x, minY + y, b): " + r.getSample(minX + x, minY + y, b) );
							totals[n][b]= totals[n][b] + r.getSample(minX + x, minY + y, b) ;
						}
					}
					//totals[n][b] = totals[n][b] / totalsMax[b];	
				}//if isBinary = false
				//System.out.println("IqmFracPyramidOperator: n: " + n + "totals[n][b]: "+ totals[n][b] );
				//totals[n][b] = totals[n][b]/totalsMax[b]; //values are already normalized
				//System.out.println("IqmFracPyramidOperator: totals[n][b]  "+ totals[n][b] );
			}//b bands		
		}
		//Log
		double[][] lnTotals  = new double[number][numBands];
		double[][] lnNumbers = new double[number][numBands];
		for (int n = 0; n < number; n++){
			for(int b = 0; b < histogram.getNumBands(); b++){
				if (totals[n][b] == 0) totals[n][b] = Float.MIN_VALUE;  //to ensure that log is valid
				lnTotals[n][b]  = Math.log(totals[n][b]);
				lnNumbers[n][b] = Math.log(eps[n][b]);
				System.out.println("IqmFracPyramidOperator: n:" + n + " eps:  "+ eps[n][b] );
				System.out.println("IqmFracPyramidOperator: n:" + n + " totals[n][b]: "+ totals[n][b] );
			}
		}
		
		int numColumns = model.getColumnCount();
		model.addColumn("Dp");
		model.addColumn("StdDev");
		model.addColumn("r2");
		model.addColumn("adjustet r2");
		for (int n = 0; n < number; n++) model.addColumn("DataX_"+(n+1));
		for (int n = 0; n < number; n++) model.addColumn("DataY_"+(n+1));

		boolean isLineVisible = false;
		for(int b = 0; b < histogram.getNumBands(); b++){ //mehrere Bands
			//Plot  //nur ein Band!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			//double[] dataX = new double[number];
			//double[] dataY = new double[number];
			Vector<Double> dataX = new Vector<Double>();
			Vector<Double> dataY = new Vector<Double>();
			for (int n = 0; n < number; n++){	
					dataY.add(lnTotals[n][b]);
					dataX.add(lnNumbers[n][b]);
			}
//			System.out.println("IqmFracPyramidOperator: dataY: "+ dataY);
//			System.out.println("IqmFracPyramidOperator: dataX: "+ dataX);
			if (optShowPlot){
				IqmTools.displayRegressionPlotXY(dataX, dataY, isLineVisible,
												 imgName+"_Band"+b, "Pyramid Dimension", "ln(1/2^n)", "ln(Count)",
												 regStart, regEnd, optDeleteExistingPlot);
			}
			double[] p = IqmTools.getLinearRegression(IqmTools.reverseVector(dataX), IqmTools.reverseVector(dataY), regStart, regEnd);
			//double[] p = IqmTools.getLinearRegression(dataX, dataY, regStart, regEnd);
			
			double dim = 0.0;
			if (isBinary) {
				dim = p[1];
			}else{
				//dim = 3-(p[1]/2);
				//dim = 2.0+p[1];
				dim = p[1]; //?
			}		
			model.setValueAt(dim,   b, numColumns);
			model.setValueAt(p[3],  b, numColumns+1);	
			model.setValueAt(p[4],  b, numColumns+2);	
			model.setValueAt(p[5],  b, numColumns+3);	
			//set regression data
			for (int n = 0; n < number; n++){	
				//model.setValueAt(lnEps[n][b],  b, numColumns+3+(n+1));
				model.setValueAt(eps[n][b],  b, numColumns+3+(n+1));
			}
			for (int n = 0; n < number; n++){	
				//model.setValueAt(lnTotals[n][b],  b, numColumns+3+number+(n+1));
				model.setValueAt(totals[n][b],  b, numColumns+3+ number+(n+1));
			}
		}	
		model.addTableModelListener(jTable);
		model.fireTableStructureChanged(); //this is mandatory because it updates the table
		jTable.getColumnModel().getColumn(2).setPreferredWidth(30);  //Band					  
		jTable.getColumnModel().getColumn(3).setPreferredWidth(30);  //RegStart				  
		jTable.getColumnModel().getColumn(4).setPreferredWidth(30);  //RegEnd					  	  
		jTable.getColumnModel().getColumn(5).setPreferredWidth(30);  //Interp					  	  
        return jTable;
	    //return ob;   //Optional wird ein Bild ausgegeben
  }
}
